Make your satellite image pop!¶
In [1]:
import re
import numpy as np
import matplotlib.pyplot as plt
import boto3
import rasterio
from rasterio.plot import show, show_hist
from skimage.exposure import equalize_adapthist
In [3]:
session = boto3.Session(profile_name='default')
In [4]:
scenes = [
'LC08_L1TP_205021_20180625_20180704_01_T1',
'LC08_L1TP_205021_20180727_20180731_01_T1'
]
In [5]:
def get_aws_band(scene, band, profile=False):
path, row = re.match(r'LC08_L1TP_(\d{3})(\d{3})', scene).groups()
prefix = f'c1/L8/{path}/{row}/{scene}'
with rasterio.Env(session=session):
with rasterio.open(
f's3://landsat-pds/{prefix}/{scene}_B{band}.TIF') as src:
if profile:
return src.profile
return src.read(out_shape=(src.height//10, src.width//10))
In [6]:
stacked_scenes = []
for scene in scenes:
red = get_aws_band(scene, 4)[0]
green = get_aws_band(scene, 3)[0]
blue = get_aws_band(scene, 2)[0]
stacked_scenes.append(np.stack([red, green, blue]))
Let's plot our downloaded image¶
In [11]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(stacked_scenes[0], title="Scene 1", ax=axes[0])
show(stacked_scenes[1], title="Scene 2", ax=axes[1])
Out[11]:
In [12]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,65535)
axes[1].set_xlim(0,65535)
axes[0].set_ylim(0,140000)
axes[1].set_ylim(0,140000)
show_hist(
stacked_scenes[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
stacked_scenes[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
In [13]:
!gdalinfo data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif
In [14]:
!gdalinfo data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif | grep "Band"
Ok, let's just make it 8bit, that'll work, right...¶
In [15]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(stacked_scenes[0].astype(np.uint8), title="Scene 1", ax=axes[0])
show(stacked_scenes[1].astype(np.uint8), title="Scene 2", ax=axes[1])
Out[15]:
In [16]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,50000)
axes[1].set_ylim(0,50000)
show_hist(
stacked_scenes[0].astype(np.uint8), bins=50, lw=0.0, stacked=False,
alpha=0.3, histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
stacked_scenes[1].astype(np.uint8), bins=50, lw=0.0, stacked=False,
alpha=0.3, histtype='stepfilled', title="Histogram", ax=axes[1])
That's pretty psychedelic, but not helpful¶
In [17]:
def linear_rescale(image, in_range=[0, 65535], out_range=[1, 255]):
"""
Linear rescaling
"""
image = np.clip(image, in_range[0], in_range[1]) - in_range[0]
image = image / float(in_range[1] - in_range[0])
return (
image * (out_range[1] - out_range[0]) + out_range[0]
).astype(np.uint8)
In [18]:
linear_rescaled_images = [linear_rescale(stack) for stack in stacked_scenes]
In [19]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(linear_rescaled_images[1], title="Scene 2", ax=axes[1])
Out[19]:
In [20]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
linear_rescaled_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
linear_rescaled_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
In [21]:
linear_rescaled_images = [
linear_rescale(stack, in_range=[0, 14000]) for stack in stacked_scenes
]
In [22]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(linear_rescaled_images[1], title="Scene 2", ax=axes[1])
Out[22]:
In [23]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
linear_rescaled_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
linear_rescaled_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
That's not very pretty¶
In [24]:
def percentile_stretch(image, lower_perc=2, upper_perc=98):
mask = np.all(image != 0, axis=0).astype(np.uint8) * 255
image_out = np.zeros_like(image).astype(np.uint8)
for band in range(image.shape[0]):
input_range = np.percentile(
image[band][mask > 0],
(lower_perc, upper_perc)
).tolist()
image_out[band] = np.where(
mask,
linear_rescale(image[band],
in_range=input_range,
out_range=[0, 255]),
0
)
return image_out
In [25]:
percentile_stretched_images = [
percentile_stretch(stack) for stack in stacked_scenes
]
In [26]:
fig, axes = plt.subplots(1, 2, figsize=(30,15))
show(percentile_stretched_images[0], title="Scene 1", ax=axes[0])
show(percentile_stretched_images[1], title="Scene 2", ax=axes[1])
Out[26]:
In [27]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,80000)
axes[1].set_ylim(0,80000)
show_hist(
percentile_stretched_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
percentile_stretched_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
Give me more detail¶
In [28]:
def clahe_equalize(image, kernel_size=2000, clip_limit=0.03):
image_out = np.zeros_like(image).astype(np.float64)
for band in range(image.shape[0]):
image_out[band] = equalize_adapthist(
image[band], kernel_size=kernel_size, clip_limit=clip_limit
)
return image_out
In [29]:
raster_equalized_images = [
clahe_equalize(stack, clip_limit=0.05) for stack in stacked_scenes
]
In [30]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(
(raster_equalized_images[0]* 255).astype(np.uint8),
title="Scene 1", ax=axes[0]
)
show(
(raster_equalized_images[1]* 255).astype(np.uint8),
title="Scene 2", ax=axes[1]
)
Out[30]:
In [31]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,80000)
axes[1].set_ylim(0,80000)
show_hist(
(raster_equalized_images[0]* 255).astype(np.uint8),
bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled',
title="Histogram", ax=axes[0]
)
show_hist(
(raster_equalized_images[1]* 255).astype(np.uint8),
bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled',
title="Histogram", ax=axes[1]
)
We could also remove the atmosphere...kind of...¶
In [32]:
def find_dark_object_value(arr):
"""Find the value of the dark object
ie the first dark value that is not an outlier
"""
preval = None
step = arr.max() / 255.0
for val in np.unique(arr)[:100]:
if val == 0:
continue
if preval is not None and (val - preval) < step:
break
else:
preval = val
return preval
def dos(image):
red = image[0]
green = image[1]
blue = image[2]
# get the dark object for each band
dark_r = find_dark_object_value(red)
dark_g = find_dark_object_value(green)
dark_b = find_dark_object_value(blue)
corr_red = red - dark_r.astype(np.float)
corr_green = green - dark_g.astype(np.float)
corr_blue = blue - dark_b.astype(np.float)
# correct negative values to zero
corr_red[corr_red < 0] = 0
corr_green[corr_green < 0] = 0
corr_blue[corr_blue < 0] = 0
return np.stack([corr_red, corr_green, corr_blue]).astype(np.uint8)
In [33]:
prepped_scenes = [linear_rescale(stack, [0, 14000]) for stack in stacked_scenes]
dos_scenes = [dos(prepped_scene) for prepped_scene in prepped_scenes]
In [43]:
fig, axes = plt.subplots(2,2, figsize=(10,10))
show(linear_rescaled_images[0], title="Scene 1", ax=axes[0, 0])
show(dos_scenes[0], title="Scene 2", ax=axes[0, 1])
show(linear_rescaled_images[1], title="Scene 1", ax=axes[1, 0])
show(dos_scenes[1], title="Scene 2", ax=axes[1, 1])
Out[43]:
In [35]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
dos_scenes[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
dos_scenes[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
Export for display¶
In [36]:
def write_rendered_rgb(rgb, profile, outname_prefix):
"""
Write a rendered RGB to GeoTIFF.
"""
profile.update(dtype=rasterio.uint8, count=3, photometric='RGB')
outname = f'data/{outname_prefix}_rendered.tif'
with rasterio.open(outname, 'w', **profile) as dst:
for k, arr in enumerate(rgb):
dst.write_band(k+1, arr)
In [37]:
profile = get_aws_band(scenes[1], 4, True)
equalized_image_8bit = (raster_equalized_images[1] * 255).astype(np.uint8)
write_rendered_rgb(equalized_image_8bit, profile, "scene2")
In [38]:
!gdalinfo data/scene2_rendered.tif
In [39]:
!gdalinfo data/scene2_rendered.tif | grep Band